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We report results of molecular dynamics simulations of amorphous ice for pressures up to 22.5 
kbar. The high-density amorphous ice (HDA) as prepared by pressure-induced amorphization of Ih 
ice at T — 80 K is annealed to T = 170 K at various pressures to allow for relaxation. Upon increase 
of pressure, relaxed amorphous ice undergoes a pronounced change of structure, ranging from the 
low-density amorphous ice (LDA) at p = 0, through a continuum of HDA states to the limiting 
very high-density amorphous ice (VHDA) regime above 10 kbar. The main part of the overall 
structural change takes place within the HDA megabasin, which includes a variety of structures 
with quite different local and medium-range order as well as network topology and spans a broad 
range of densities. The VHDA represents the limit to densification by adapting the hydrogen-bonded 
network topology, without creating interpenetrating networks. The connection between structure 
and metastability of various forms upon decompression and heating is studied and discussed. We 
also discuss the analogy with amorphous and crystalline silica. Finally, some conclusions concerning 
the relation between amorphous ice and supercooled water are drawn. 
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I. INTRODUCTION 

Properties of amorphous solid water at low tempera- 
tures continue to attract the attention of both theory and 
experiment. It has been known for a long time that at 
least two distinct amorphous forms of ice exist. High- 
density amorphous ice (HDA) is prepared by compres- 
sion of ordinary Ih ice to 12 kbasi and when recovered 
at ambient pressure it has a density of ~ 1.17 g/cm 3 . 
Upon isobaric heating to 117 K, the density drops con- 
siderably and a second form is found, called low-density 
amorphous ice (LDA)i with a density ~ 0.94 g/cm 3 . The 
transition between LDA and HDA can also be induced 
by pressure^. Interest in this system is enhanced by the 
possible existence of a second critical point in water, pro- 
posed originally in Ref. 4 and later elaborated in several 
variants (RefsA^ji). According to this hypothesis, in the 
deeply supercooled region a second critical point should 
exist, below which two kinds of water, high-density liquid 
(HDL) and low-density liquid (LDL) are separated via a 
first-order phase transition. Experimentally, however, it 
has not yet been possible to access the deeply supercooled 
region and directly investigate its properties. In the lack 
of direct evidence, the existence of two different amor- 
phous forms of ice has been used as an indirect support 
for this hypothesis, assuming that HDA and LDA, appar- 
ently separated by a sharp transition, are simply glassy 
forms of HDL and LDlA 

Recently several new experimental 

results£*i2iii*i2i±2ii2ii£ii& raised new questions about 
the nature of amorphous ices as well as about their rela- 
tion to the speculated second critical point of water. A 
new amorphous form of ice has been reported—, prepared 
by heating HDA under pressure of 11 kbar to T ~ 170 
K and cooling it back to T = 80 K; when recovered 



at ambient pressure it has a density of ~ 1.25 g/cm 3 . 
It has been called very high-density amorphous ice 
(VHDA) and characterized experimentally by neutron 
diffractionii. More detailed structural measurements of 
VHDA ice were performed very recentlj*i£ and showed 
that the radial distribution function (RDF) of VHDA 
is in fact more structured than that of HDA and LDA, 
revealing the presence of an enhanced medium-range 
order. Other experimentsi2iiiiiii& focused on the HDA- 
LDA transition induced by heating at ambient or low 
pressure. In Ref»i£ it was shown that by heating HDA 
to temperatures intermediate between 80 K and 110 K 
the sample can be trapped in an apparent continuum 
of metastable structures between HDA and LDA. This 
suggests that there might be no sharp transition between 
the two forms. In Refii 3 - the kinetics of the HDA to 
LDA transition was studied; this revealed three stages of 
the transformation, consisting of the annealing of HDA, 
followed by an accelerated transition to the LDA form 
and subsequent annealing of the LDA. The experiment 
in Ref. 15 , while also finding a continuum of HDA states, 
observed also a propagation of the LDA-HDA interface, 
thus providing evidence in favor of a sharp transition 
between the two forms. Possible implications of the new 
experiments have been discussecUSiiSiiS. Several review 
papers on amorphous and supercooled water have also 
appeared recently, see Refs^jSi^. 

Simulations can complement the experiment by pro- 
viding access to shorter time scales, not easily accessed 
in experiment. At the same time they can provide de- 
tailed structural information which might not be easily 
extracted from the experimental data. The basic phe- 
nomenology of amorphous ices has been reproduced in 
earlier works&2ii24. New simulations have also been per- 
formed recentlj&^KA^- 
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Here we present the results of extensive constant- 
pressure MD simulations of amorphous ice at tempera- 
tures 80 - 170 K and pressures up to 22.5 kbar. We focus 
on the structure of the annealed form of amorphous ice 
and study its evolution under increasing pressure. We 
identify the LDA and HDA regimes, possibly separated 
by a transition. VHDA on the other hand does not ap- 
pear to be a new phase, but rather a particular high- 
pressure regime of HDA. In particular, we suggest that 
both new phenomena, the VHDA and the continuum of 
metastable HDA densities at ambient pressure, originate 
from the relationship between the density and the topol- 
ogy of the hydrogen-bonded network of the HDA phase. 
A preliminary account of this work has already been pub- 
lished in Refill. Here we present a more detailed account 
of the results as well as new simulations and new anal- 
ysis. The paper is organized as follows. In section [H] 
we present the model and details of our simulation tech- 
nique. In section II I II we present our results. First we 
describe the protocol used to prepare various forms of 
amorphous ice and then discuss in detail their properties, 
comparing with experiment as well as with other simu- 
lations. In section IIVI we draw an analogy of another 
important tetrahedrally-bonded system, amorphous sil- 
ica. Finally, in the last sectionEJwe summarize our find- 
ings and draw some conclusions concerning the relation 
between amorphous ice and supercooled water. 

II. MODEL AND TECHNIQUE 

Our tool is constant-pressure molecular dynamics 
(MD) simulation. We employed the GROMACS code 32 
using the Parrinello- Rahman constant-pressure MD 33 , 
the Berendsen thermostat^ and a time step of 2 fs. Elec- 
trostatic interactions were treated by the particle-mesh 
Ewald method— An initial proton-disordered configu- 
ration of Ih ice with 360 H2O molecules and zero dipole 
moment in an orthorhombic box was prepared using the 
Monte Carlo procedure of Ref— . 

The interaction between water molecules was described 
by the classical TIP4P model 37 . In previous simula- 
tions this was found to reproduce well the transitions 
I h - HDA and LDA - HDA, both qualitatively and 
quantitativelySiSiSi. Although it is known to system- 
atically overestimate the experimental ice densities by 
about 2 %, the TIP4P model has recently been shown to 
be capable of predicting qualitatively correctly the entire 
phase diagram of water—. 

III. RESULTS AND DISCUSSION 

A. Preparation and annealing of amorphous ice 

First we shall describe the simulation protocol that 
we applied. We started by compressing the ice Ih at 
T = 80 K, increasing the pressure in steps of 1.5 kbar. 



At 13.5 kbar a sharp transition occurs and the density 
increases by almost 30 % to 1.37 g/cm 3 (Fig|TJ. The 
sample was then further compressed at T — 80 K to 22.5 
kbar and from 15 kbar decompressed to p = 0. This par- 
ticular HDA form obtained by pressure-induced amor- 
phization of Ih ice at T = 80 K and subsequent compres- 
sion or decompression, without any thermal treatment, 
will be in the following denoted by HDA' (as-prepared 
HDA phase). During decompression the HDA' density 
gradually decreased and at p — reached the value of 
1.19 g/cm 3 , close to the HDA experimental value of 1.17 
g/cm 3 i The radial distribution function (RDF) of HDA' 
at p = is shown in Fig|5| it has a broad second peak 
between 3.3 and 4.6 A, very similar to that of HDA at 
p = C 10 . Inspired by the experiments^ that led to the 
discovery of the VHDA we decided to anneal the HDA' 
phase at each intermediate pressure between 22.5 kbar 
and zero (a similar treatment was applied in Ref£ at 
p = 8.4, 11 and 19 kbar) in order to search for possible 
new structures. Annealing was performed by heating up 
to T = 170 K and subsequent cooling down to 80 K; the 
temperature was always changed in steps of 10 K. The 
phase obtained in this way will be called relaxed phase 
(RP). We note here that in Ref— a slow molecular dif- 
fusion was observed in the MD simulation of the TIP4P 
model at T — 173 K; therefore at 170 K it should be pos- 
sible to equilibrate the system within accessible simula- 
tion times. While in experiments a HDA' sample heated 
at an arbitrary pressure might recrystallize&iS^i, in the 
time scale of a simulation this is not likely to happen. 
We are therefore restricted to exploring the (metastable) 
disordered structures. In order to check the metastability 
of the RP phase upon decompression, the RP prepared 
at each pressure was finally decompressed at T — 80 K 
down to p = 0, decreasing the pressure in steps of 1.5 
kbar. 
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FIG. 1: Density vs. pressure dependence at T = 
80 K for the various amorphous phases during compres- 
sion/decompression. The triangles point in the direction of 
pressure change, the lines are just a guide for the eye. 
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FIG. 2: Oxygen - oxygen radial distribution function of var- 
ious amorphous phases at T — 80 K and p = — 22.5 kbar. 
The curves in the upper part of the figure have been shifted 
by 8 for clarity. 



After each change of pressure or temperature the sys- 
tem was equilibrated for 5 ns and observable quantities 
were averaged over 0.5 ns. An additional equilibration 
for 25 - 50 ns was performed during the annealing at the 
highest temperature of T = 170 K. We stress here that 
the long equilibration times are indeed necessary in order 
to allow the system to undergo structural changes; e.g at 
pressure p — 0.75 kbar the equilibration of the system at 
T = 170 K takes about 20 ns. A comprehensive summary 
of the density vs. pressure dependence of Ih, HDA', RP 
and decompressed RP phases at T = 80 K is presented in 
Fig^and will be discussed in detail in the next section. 

In Fig[3]we show the relaxation of the density during 
annealing of the HDA' ice at p = 16.5 kbar to T = 170 K. 
It can be seen that upon annealing at 90 and 100 K the 
density grows while between 110 and 130 K the sample 
undergoes a thermal expansion. At 140 K and 150 K 
the density grows further but no substantial change is 
observed above T = 150 K. 

We have verified that the enthalpy of the relaxed phase 
is lower than that of HDA' ice at any pressure (FigQJ. 
Assuming that the entropy contribution to the Gibbs po- 
tential can be neglected at T = 80 K (and entropy differ- 
ences between amorphous phases are anyway expected to 
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FIG. 3: Time dependence of the density during annealing 
of the HDA' ice at p = 16.5 kbar to T = 170 K. Points 
where the temperature is increased by 10 K are marked by 
vertical dashed lines. The temperature (in K) during each 
time interval is shown in the top part of the figure. 



be small^i) this demonstrates that upon annealing at any 
pressure HDA' ice irreversibly relaxes to a state with a 
lower free energy. This a posteriori justifies the necessity 
of annealing in order to reach a metastable equilibrium 
corresponding at each pressure to a well-defined ther- 
modynamic phase. We note that the lowest difference 
between the enthalpy of the HDA' phase and that of the 
RP is found at p = 4.5 kbar, suggesting that at the lat- 
ter pressure HDA' is closest to its corresponding relaxed 
amorphous form (see also next subsection). 
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FIG. 4: Decrease of enthalpy upon annealing: difference be- 
tween the enthalpy of the RP phase and that of the HDA' 
phase at T=80 K. 



B. Evolution of the RP with increasing pressure 

In this section we analyze the properties of the RP and 
their dependence on pressure, with the focus on struc- 
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ture. The density vs. pressure dependence of the RP 
(Fig. can be considered as the equation of state of 
amorphous ice. We note first that the HDA' and RP 
curves cross at about 7 kbar; for lower pressure anneal- 
ing results in lower density while for higher pressures the 
system densities. 

At p = the density after annealing reaches a value 
of 0.97 g/cm 3 , which agrees well with the experimental 
LDA value of 0.94 g/cm 3 . The remarkable feature of 
the RP curve is the narrow region between 1.5 and 2.25 
kbar where the density increases by about 9 %, from 1.04 
g/cm 3 to 1.13 g/cm 3 . Upon further increase of pressure 
the density grows fast and reaches at p = 3 kbar the value 
of 1.19 g/cm 3 . Beyond that point the density growth 
slows down progressively and at the highest pressure of 
22.5 kbar the density reaches a value of 1.49 g/cm 3 . 

The 0-0 RDF's of RP at different pressures are shown 
in Fig|3 We also calculated the O-H RDF (not shown) 
for RP at p = 0,2.25,6 and 15 kbar. Integrating be- 
tween 1.5 and 2.25 A we found at all pressures a coordi- 
nation number of 2, indicating a fully hydrogen-bonded 
network. In order to characterize the evolution of the net- 
work topology we calculated the ring statistics for RP at 
all pressures. This reveals information on medium-range 
order that otherwise might not be easily extracted from 
the RDF^Sii 3 .. We applied the ring counting algorithm 44 
using the criterion from Ref4£ to identify the hydrogen 
bonds (we used for the 0-0 distance a cutoff parameter 
of r cut — 3.05 A and a tolerance of A = 0.2 A). Our 
aim is to draw qualitative conclusions, so we did not try 
to bring down the statistical error by repeating several 
times the quench from 170 K to 80 K and averaging over 
the resulting structures. We now discuss the evolution 
of the RP at increasing pressure in terms of RDF and 
network ring statistics (Fig[5J| and show that there are 3 
distinct regimes. 

The RDF of the RP at p = (Fig© exhibits at r = 3.1 
A a very deep minimum between the first and second 
shell and a well-defined second shell peak at r = 4.4 
A, very similar to that found experimentally for LDA 
in Ref4&. The phase thus coincides with the LDA as 
expected. The same is true at 0.75 kbar at a density 
of 0.99 g/cm 3 , where the RDF within the second peak 
is practically indistinguishable from the one at p = 0, 
and only beyond 5 A can small differences be seen. At 
1.5 kbar the density increases to 1.04 g/cm 3 while the 
RDF becomes slightly different from that of LDA at p = 
0; the very deep minimum between the first and second 
shell is still present. At all these pressures the network 
is dominated by 6-membered rings. 

The properties of the RP at p — 2.25 kbar are rather 
different, correlating with the sharp increase of density. 
The second shell peak of the RDF drops and shifts to 
lower r and at the same time the RDF grows substan- 
tially in the region around r = 3.3 A, revealing the pres- 
ence of interstitial molecules 4 - . A dramatic change is 
seen in the ring statistics: the number of 6-membered 
rings now starts to drop and at the same time the num- 
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FIG. 5: Number of n-membered network rings in the RP as 
a function of pressure at T = 80 K and p = — 22.5 kbar. 



ber of 8-membered rings grows fast. This behavior is 
compatible with a transition from LDA to HDA occur- 
ring between 1.5 and 2.25 kbar. Around p = 4.5 kbar the 
RDF develops a broad second peak between 3.2 and 4.4 
A, similar to that of HDA' at p = (Fig©. This agrees 
with the observation based on the enthalpy difference and 
shows that the HDA' is indeed closest to the RP at this 
pressure. Approaching p ~ 10 kbar the ring statistics 
are definitely dominated by 8 and 9-membered rings; the 
network has thus undergone a substantial reconstruction 
(see also section ITV|) . For p > 10 kbar, the ring statis- 
tics almost stabilize, revealing that the reconstruction of 
the network is practically completed; at the same time 
the density growth slows down further and the compress- 
ibility approaches that of ice Ih- This indicates that the 
increase of density due to the more efficient packing of 
the molecules has at p ~ 10 kbar reached its limit at 
the value of pu m ^1.39 g/cm 3 , and further compression 
proceeds mainly by elastic compression. In this limiting 
regime the RDF develops a pronounced second peak at 
r = 3.2 A (FigEJl while the original second shell peak at 
r = A A A disappears completely. In the subsection IIII El 
we will identify this regime with the VHDA form. 

The above analysis provides quite a clear picture of 
the evolution of the network topology when going from 
LDA to VHDA. Concerning this point, rather contra- 
dictory opinions have been presented in the literature, 
based on the indirect information provided by the de- 
tailed analysis of radial and spatial distribution functions 
obtained from diffraction experiments and the empirical 
potential structure refinement procedures^. First we note 
that there is no sign of any discontinuity upon evolution 
of the HDA phase into the limiting VHDA regime. In 
Refiii it was speculated that there is no need to postu- 
late any significant reorganization of the network struc- 
ture in moving between HDA and VHDA and that they 
appear topologically isomorphous. Our results, however, 
show clearly that the evolution of HDA between 2.25 kbar 
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and 10 kbar must inevitably involve a substantial recon- 
struction of the network. On the other hand, in Ref. 10 
it was suggested that HDA' under pressure shows some 
characteristics of interpenetrating networks, similar to 
those of high-pressure crystalline ices VII and VIII. We 
checked this feature in the RP up to the highest pres- 
sure investigated, using the following algorithm. Start- 
ing from each molecule we considered a sphere with ra- 
dius r cut = 5 A and tested whether the molecules within 
the sphere were connected to the central one by a path 
containing no more than 4 hydrogen bonds. We found 
that practically all molecules fulfilled this criterion; this 
is clearly incompatible with the presence of interpene- 
trating networks^. The VHDA structure can therefore 
be considered as the upper limit to efficient amorphous 
packing of the molecules without creating interpenetrat- 
ing networks, as suggested in Ref?^. Still, it is an in- 
teresting question whether amorphous ice with interpen- 
etrating networks can exist. Very recently, a study of 
VHDA was performed^ in the region of densities rang- 
ing up to 1.9 g/cm 3 , thus much higher than those stud- 
ied here. Based on the analysis of bond angle distri- 
butions, they showed that VHDA at very high densities 
approaches a disordered close-packed structure, with lo- 
cal order more similar to the fcc/hcp than to the bcc 
crystal. From this they concluded that VHDA does not 
represent a disordered version of ice VII and therefore 
does not have interpenetrating networks. The algorithm 
applied here in fact provides direct evidence. We have 
also generated a sample of the RP at p — 42 kbar, with 
density p — 1.61 g/cm 3 , and found that the ring statis- 
tics were practically equal to the average of the RP in 
the interval 10 - 22.5 kbar. This shows that the network 
topology stabilizes already for densities p > pn m ^1.39 
g/cm 3 , although the local order converges only at much 
higher densities. 



C. Analysis of the shape of the rings 

In order to obtain a deeper insight into the structural 
response to pressure of the RP in the three regimes (LDA, 
HDA, VHDA), we also performed a more detailed analy- 
sis of the network structure, separating the contributions 
coming from different rings. 

The shape of the rings has been characterized through 
the eigenvalues of the inertia tensor. Given the three 
eigenvalues I\, I2 and 23, sorted by increasing magnitude, 
we define an elongation parameter east= (1% — 1\) / 13 
and the asphericity (a) as the root mean square devia- 
tion of Ii normalized to I3. According to the definitions, 
e can vary from 0, for a circular ring or a sphere, to 1 
for a linear arrangement, while a is zero for a sphere, 
0.236 for a circular ring and 0.58 for a linear arrange- 
ment. The distribution of a for all the rings considered 
is peaked around ~ 0.2 at p = 0. Increasing the pres- 
sure increases the spreads of the distribution and shifts 
the peak toward higher a. The tail of the distribution 



never extends below 0.15, which means that the rings are 
mainly planar structures. 

The probability distribution of the elongation param- 
eter provides a clear indication of the evolution of the 
shape of different rings as a function of the pressure. 
In HDA at p = 2.25 kbar (FigEJi) the P(e) shows that 
the water molecules in five and six-membered rings ar- 
range themselves in circular rings. The larger the ring 
the broader the distribution, indicating that larger rings 
can arrange into elongated structures without paying too 
much in terms of strain energy. P(e) at different pres- 
sures for six and nine-membered rings are shown in FigO 
panels (b) and (c), respectively. Six and nine-membered 
rings, the quantity of which is most affected by pressure- 
induced phase transitions, are shown as representative of 
the behavior of small and large rings under pressure, re- 
spectively. At low pressures (up to 2.25 kbar) the shape 
of the small rings (FigHJJ)) remains unchanged but when 
the transition to HDA occurs their number rapidly de- 
creases, and elongated large eight and nine-membered 
rings are formed. In fact, already at p = 1.5 kbar the 
main peak of the -P(e) of nine-membered rings (FiglU) 
shifts from 0.1 to 0.48. The broad P(e)'s of large rings at 
higher pressures show that such rings can easily adapt to 
any shape so as to achieve a more efficient compaction. In 
the HDA region the morphology of small rings evolves to- 
ward more elongated shapes and their amount decreases 
continuously. When the onset pressure for the VHDA 
regime is reached, the topology of the amorphous net- 
work stops changing and also the P (e)'s stabilize. The 
residual small rings in VHDA are arranged into elongated 
structures (the P(e) is peaked at 0.4), whereas the P(e) 
of the nine-membered rings has no sharp peaks. At this 
regime the response to further compression consists in 
the deformation of short-range structural features, such 
as the bond angle distributions, as no more important 
rebonding is induced by increasing the pressure up to 42 
kbar. 

We define a ring-restricted radial distribution function 
n-rRDF(r) as the probability of finding two atoms at 
a distance r within the same n-membered ring. This 
quantity allows us both to recognize the separate contri- 
bution of different rings to the g(r) and to characterize 
the response of different rings to compression. In Fig[7| 
the oxygen-oxygen n-rRDF(r) at different pressures are 
shown for six and nine-membered rings. The position 
of the first peak of both rRDF's is unaffected by com- 
pression up to 15 kbar and even a compression to 42 
kbar induces a shift as small as 0.04 A. Low pressures 
up to 2.25 kbar do not affect significantly the 6-rRDF, 
whereas in the same range of pressures nine-membered 
rings already provide a sizable contribution to the inter- 
stitial region between the first and the second peak of 
the g(r). In the region of stability of the HDA (6 kbar) 
six-membered rings are strained and contribute to the 
interstitial region of the g{r) through a broadening and a 
shifting to the left of the second peak. On the other hand 
the second peak of the 9-rRDF(r) spreads partly in the 
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FIG. 7: Ring-restricted radial distribution function at several 
pressures for six (a) and nine-membered (b) rings. 



interstitial region and gives rise to a shallow peak at 4.9 
A. Such a peak becomes more pronounced as the pressure 
increases and is a feature of VHDA 17 . It is worth not- 
ing that even at 42 kbar the six-membered rings do not 
provide a sharp contribution to this feature. We note 
that the rRDFs fully account for the interstitial peak 
in the g(r) of VHDA, showing that this peak originates 
from contributions within the same ring. This consti- 
tutes independent evidence that there is no formation of 
interpenetrating networks. 



FIG. 6: (a) Probability distribution of the elongation param- 
eter of the rings in HDA at 2.25 kbar. (b) and (c) P(e) at sev- 
eral pressures from LDA to VHDA for six and nine-membered 
rings, respectively. We note that all the curves are normal- 
ized to one and therefore do not reflect the change of the total 
number of rings with a given size. 



D. Transition LDA-HDA 

The question of whether there is indeed a transition 
between LDA and HDA and, if so, what its character is, 
is of great importance. From the ring statistics shown 
in Fig|S] it can be seen that the response of the num- 
ber of rings to the applied pressure is rather different in 
the regions below 1.5 and above 2.25 kbar. Because of 
the limited accuracy of the ring statistics calculated in 
the frozen states, we can conclude that this behavior is 
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compatible with the existence of a phase transition as 
required by the second critical point scenario^. 

In order to determine whether the density of the RP 
changes discontinuously between 1.5 and 2.25 kbar, we 
performed at T — 170 K also a ~ 30 ns simulation at an 
intermediate pressure of 1.875 kbar. We did not observe 
any mctastability or hysteresis effects but instead rather 
large and slow density fluctuations. Since our system 
is relatively small, such behavior is compatible with a 
weak first order transition, occurring just below the crit- 
ical temperature, where the system oscillates over a low 
barrier between two states. In principle it is, however, 
also possible that the change of density with pressure is 
genuinely continuous, with a highly compressible region 
between 1.5 and 2.25 kbar. In order to shed more light 
on this important issue, it would be necessary to perform 
simulations with larger systems, including several system 
sizes, and apply standard finite-size scaling techniques 49 . 
It might also be useful to perform free energy calcu- 
lations, e.g. umbrella sampling, with a suitable order 
parameter, similar to that performed in RefiSi. Both 
techniques would, however, require considerable CPU re- 
sources to achieve a good equilibration and sampling, 
since in the interesting region the free energy surface is 
very rough, resulting in very long autocorrelation times. 



E. Decompression of RP to p = 

Most experimental data on amorphous ices have been 
gathered on samples decompressed to ambient pressure. 
To our knowledge, there are no experimental data for the 
RP under pressure to which we could directly compare 
our results of the subsection IIII Bl We discuss now the 
interesting behavior of the RP upon decompression. 
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FIG. 8: Density of the RP after decompression at T — 80 K 
from pressure p. 

The density of the LDA phase from 0.75 kbar relaxes at 
p = to p = 0.98 g/cm 3 , very close to the p = 0.97 g/cm 3 
of LDA prepared at p = 0. The LDA phase from 1.5 kbar 



on the other hand reaches at p = a somewhat higher 
value of p = 1.025 g/cm 3 , suggesting the possibility that 
at T = 80 K and p = there might not be a unique struc- 
ture of the LDA (in agreement with RefsiS&fi 3 -) and this 
phase can actually span a narrow density interval. In the 
pressure interval p = 2.25 — 10 kbar, all decompression 
curves are roughly parallel (Fig[Q and the faster growing 
RP density results upon decompression in an increasing 
density at p = 0, spanning the interval from 1.10 to 1.26 
g/cm 3 . The picture changes remarkably for p > 10 kbar. 
Here, the slope of the RP curve becomes close to that 
of the decompression curves which lie close to each other 
and initially almost follow the RP curve. Decompression 
from almost all pressures results at p — in a narrow 
interval of densities around pvhda ~ 1-29 g/cm 3 , corre- 
sponding to the decompression from the limiting density 
piim- For convenience, in Fig[S]we show the dependence 
of the final p = density on the original pressure p at 
which the RP was prepared, where the saturation can 
be clearly seen. This agrees very well with the obser- 
vation in Ref. 9 where the samples annealed at 11 and 
19 kbar reached upon decompression the same VHDA 
density of 1.25 g/cm 3 ; in fact, this was the reason why 
VHDA was originally suspected to represent a new ther- 
modynamic phase. The RDF of RP decompressed from 
15 kbar (Fig(2J is clearly similar to that of VHDA recov- 
ered at p = in experiment 11 , showing the presence of 
the distinct peak at r — 3.37 A, very close to the first shell 
peak. This allows us to identify this form as VHDA. The 
spectrum of metastable states at p = (FigJSJ is thus 
compatible with the existence of a narrow LDA region 
and a broad continuum of metastable HDA states with a 
density below that of VHDA (as found experimentally in 
Refii^). We note that the density spectrum of the HDA 
states extends both above and below that of the HDA', 
and there is no reason to consider HDA' as a particu- 
lar state representative of the HDA phase. While the 
LDA states might be separated from the HDA ones by a 
gap, it must certainly be much smaller than the density 
difference between LDA and HDA' at p = 0. 

We now suggest an explanation for the existence of an 
upper limit pvhda to the density of metastable HDA at 
T = 80 K and p = 0. In the HDA phase with p < pn m , 
the system with increasing pressure achieves a better 
packing of molecules due to network reconstruction. This 
must involve the breaking and remaking of bonds, which 
at any density requires crossing a free energy barrier, and 
is only possible due to the annealing at higher temper- 
ature, in our case 170 K. During cooling to 80 K, the 
network topology becomes frozen. Upon subsequent de- 
compression at 80 K, the barriers cannot be crossed and 
therefore the system cannot relax to a lower density via 
reconstruction of the network. This is illustrated in Fig[5] 
where the evolution of the ring statistics upon decompres- 
sion of the RP from p — 13.5 kbar is shown; it can be seen 
that no pronounced change in the network topology oc- 
curs. The largest change is seen in the number of 9 and 
8-membered rings, which decrease by about 20 % and 
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9%, respectively, while the number of other rings stays 
practically constant. The decompression thus proceeds 
dominantly via relaxation of the elastic compression and 
that is why the HDA states with p < pu m , which have 
a variety of different topologies, relax elastically to dif- 
ferent states, spanning a range in density. On the other 
hand, all HDA states with p > pu m , which have almost 
the same network topology, relax upon decompression to 
the same density pvhda- This accounts for the behav- 
ior observed in Ref£, with no need to postulate VHDA 
to be a new phase, and is also consistent with the fact 
that we did not observe any discontinuity during the evo- 
lution from HDA to VHDA. The origin of the behavior 
of VHDA upon decompression therefore appears to be 
kinetic rather than thermodynamic. 
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FIG. 9: Ring statistics during the decompression of RP at 
T = 80 K from p = 13.5 kbar to zero. 

The two points from the highest pressures of 21 and 
22.5 kbar reach upon decompression a lower density, 1.21 
and 1.23 g/cm 3 , respectively. We believe that this is 
related to the fact that our decompression is performed 
too fast, resulting in excessive elastic energy at p = 0. 
This may, in turn, allow some barriers to be crossed and 
enable a transition to a density p < pvhda- 

We comment now on the experimentsi2ii4, where 
HDA' was heated to intermediate temperatures between 
80 and 120 K and at each temperature annealed for sev- 
eral hours. On this time scale every increase of tempera- 
ture resulted in an initial slow drop of the density which 
afterwards gradually stabilized; an additional drop of the 
density could be observed only by further increase of tem- 
perature. This behavior clearly points to the fact that 
as the HDA phase approaches its low density limit, the 
barriers increase and can be overcome only at a higher 
temperature. Our picture of the structural evolution of 
the HDA phase is compatible with these experimental 
findings. It is plausible to assume that the height of the 
barriers is correlated to the amount of network recon- 
struction necessary to change the volume, and is there- 
fore related to the pressure derivatives of the number of 
As shown in sectionllll Bl the reconstruction of the 



rings. 



RP upon increasing pressure is most dramatic at the low 
density limit of the HDA spectrum, and with increasing 
density becomes gradually less pronounced, until it prac- 
tically vanishes for p > pu m . We explicitly checked the 
degree of metastability of HDA forms with different den- 
sities (RP decompressed from several different pressures) 
by heating at p = 0; the results are shown in FigllOl 
The temperature was increased in steps of 10 K, spend- 
ing 5.5 ns at each temperature, and we note that also 
here it would be useful to perform the heating several 
times and average in order to improve the statistical er- 
ror. Nevertheless, it can be clearly seen that the most 
stable structure under heating is actually the lowest den- 
sity HDA (RP decompressed from 2.25 kbar) which un- 
dergoes only a small drop in density up to 130 K. With in- 
creasing initial density the samples start to relax at lower 
temperature, which confirms the above relation between 
the network topology and barriers. It is also seen that 
HDA' is the least stable of all the samples, as may be ex- 
pected for an insufficiently equilibrated phase possessing 
an excess free energy. We stress that at T = 170 K all 
curves reach practically the same density, lower than 0.99 
g/cm 3 , which also agrees with the experimental finding 
in Refsi^jii that VHDA upon heating converts to LDA. 
This is different from what found in Rcf. 28 where it was 
argued that VHDA upon heating converts to a different 
form of LDA, with a higher density of about 1.04 g/cm 3 . 
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O-ORP from 4.5 kbar 

O-O RP from 10.5 kbar (VHDA) 

O-O RP from 15 kbar (VHDA) 
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FIG. 10: Density as a function of temperature during heating 
of various decompressed RP phases as well as HDA' at p — 0. 



In RefS various HDA forms were prepared by an al- 
ternative procedure, consisting of the pressure-induced 
amorphization of cubic ice at different temperatures 
ranging from 50 to 300 K. It was concluded that the 
HDA' produced by pressure-induced amorphization at 
liquid nitrogen temperature and below represents a lim- 
iting form of HDA and thus the phase spans a density 
interval between HDA' and VHDA. This procedure does 
not cover the part of the HDA spectrum that has a den- 
sity below that of HDA' and can exist at p = and 
T = 80 K in metastable form and therefore we do not 
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consider the imperfectly equilibrated HDA' to be a lim- 
iting form of HDA. 

In Ref»ii the interstitial occupancy in HDA' and 
VHDA at ambient pressure was calculated by integrating 
the 0-0 RDF between 2.3 and 3.3 A. The values of 5.0 
and 5.8 were found, respectively, corresponding to 50 % 
or almost 100 % occupancy of the "lynch pin location" . It 
was speculated that due to some unknown specific mech- 
anism only these values and no intermediate ones can be 
made stably at ambient pressure. We calculated the same 
occupancy at T — 80 K and p — in our HDA' as well as 
in the RP decompressed from all pressures (Fig lll|) . In 
HDA' we found a value of 4.9 while in the HDA branch 
of decompressed RP we found an apparently continuous 
spectrum ranging from 4.3 (from 2.25 kbar) to about 6 
(from pressures above 15 kbar). This again clearly shows 
that while VHDA indeed represents a limiting structure, 
this is not the case for HDA' whose occupancy close to 
the value of 5 can be regarded as accidental. 



this is also clearly present in our RDF and in subsection 
1111 CI we have identified its origin in the contribution of 
9-membered rings (see Fig[7{b)). 
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FIG. 11: Oxygen occupation number within 3.3 A in the RP 
decompressed from pressure p at T — 80 K. 

Recently, a neutron and X-ray diffraction study of the 
VHDA structure was reported in Reftf, showing the 
presence of at least seven well-defined shells in the RDF 
of the VHDA, extending almost to the distance of ~ 20 A. 
This reveals the presence of an enhanced medium-range 
order in the VHDA. In order to check this property, we 
also prepared a bigger VHDA sample consisting of 2880 
water molecules, allowing us to calculate the RDF up to a 
distance of 20 A. We followed basically the same protocol 
as for the 360-molecule sample and annealed the HDA' 
at p = 15 kbar, but using shorter simulation times. The 
radial distribution function Doo( r ) — 47r pr(g{r) — 1) in 
the decompressed sample is shown in Fig^] Apart from 
the height of the first peak, our result agrees well with the 
experimental one (Fig. 2(b)) in Ref>i£, and also shows at 
least seven coordination shells extending beyond ~ 16 A. 
The presence of such enhanced medium-range order is 
likely to be related to the fact that the network topol- 
ogy of VHDA is dominated by large rings. In Ref>i£ the 
existence of a well-defined shell at 5 A was pointed out; 
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FIG. 12: Oxygen-oxygen radial distribution function Doo(r) 
for the VHDA prepared by annealing at p — 15 kbar and 
decompressing to zero pressure. 



IV. ANALOGY TO SILICA 

In spite of the different nature of the bonds between 
water molecules in amorphous ice and between silicon 
and oxygen atoms in amorphous silica, both systems con- 
sist of a continuous random network of corner-sharing 
tetrahedra^i and in some windows of their phase dia- 
grams display analogous phenomenologies when pressure 
is applied. Each tetrahedral unit in silica is made of a 
four-fold coordinated silicon atom in the center and four 
bridging oxygen atoms at the corners. The size distri- 
bution of the primitive rings 4 ^ 2 - is peaked at six silicon 
atoms per ring and presents a sizable amount of four 
to ten-fold ringa 42 i 44 i 52 . A high-density (HD) phase of 
amorphous silica was discovered by GrimsditchS? about 
20 years ago. It was shown that upon compression above 
8 GPa a-Si02 undergoes a permanent densification which 
amounts to about 20% when the system is released to 
ambient pressure and the transition to the HD phase is 
accompanied by irreversible structural changes, observ- 
able by Brillouin and Raman measurements^^. 

As in the case of amorphous ice, the nature of the phase 
transition is still debated, since no discontinuous volume 
change is observed in compression experimentsSiSiS 
while computer simulations do not clarify whether there 
is a kinetically hindered first order transition^ or a pres- 
sure window where there is a balance between two den- 
sification mechanisms^ 2 *^. In fact three different mi- 
croscopic mechanisms cooperate in accommodating the 
amorphous silica network in a smaller volume 4 ^^^^ 
and make the phase transition between LDA and HDA 
apparently smoother than in amorphous ice. In the low 
pressure regime (~3 GPa according to Refsi 2 -^) the 
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volume diminishes only by the buckling of the network, 
which results in a shift toward smaller values of the Si- 
O-Si bond angle distribution. In this process the tetra- 
hedral units are not deformed and no bonds are broken. 
At higher pressures the buckling mechanism is supple- 
mented by a substantial rebonding in the network, which 
mainly affects the medium-range order features: the local 
tetrahedral order is preserved but the ring-size distribu- 
tion broadens and its peak shifts to larger ringsS. In 
the response of silica to pressure the buckling and the 
rebonding mechanisms correspond to elastic and plastic 
regimes, respectively, as observed by Davila et al. in MD 
simulations^,. On the other hand, in both regimes short- 
range structural properties such as the Si-O-Si and the 
O-Si-0 angles vary continuously. In addition, coordina- 
tion defects may be formed and contribute to densifica- 
tion at even higher pressures (e.g. > 5 GPa in Ref4£). In 
ice the nature of the hydrogen bond inhibits this latter 
mechanism, as no more than four hydrogen bonds per 
water molecule can be formed: in fact the average co- 
ordination number of the RP is ~4 at all the explored 
pressures. Consequently, when the limit for the topologi- 
cal densification is reached, amorphous ice turns back to 
an elastic regime (VHDA). 

The increasing of the characteristic size of the 
rings upon densification is a general feature of both 
amorphous£2i§2i£i and crystalline^ tetrahedral networks. 
Among the tetrahedrally coordinated crystalline silica 
polymorphs, the lower density forms (cristobalite and 
tridimite) consist of six-membered rings only. In the 
denser silica crystals the average size of the rings in- 
creases accordingly to the density. For example coesite, 
which is 30% denser than a-cristobalite has an average 
ring size of ten and contains rings of size up to twelve. It 
is worth noting that unlike the HD crystalline phases of 
ice, it is impossible to form silica polymorphs with inter- 
penetrating networks, because of sterical interactions. 

V. CONCLUSIONS 

Upon increase of pressure, relaxed amorphous ice un- 
dergoes a pronounced change of structure, from LDA at 
p = to VHDA at p > 10 kbar. During this transfor- 
mation, there is possibly a discontinuity between LDA 
and HDA, although from our simulations performed on 
a relatively small system we cannot distinguish between 
a weak first-order transition and a continuous change. 
Nevertheless, we can clearly distinguish the LDA and 
HDA regimes. This identification is based on the ex- 
istence of a transition region characterized by a rapid 
change of density, on the presence of interstitial molecules 
and the behavior of the ring statistics. It is important 
to note that the main part of the overall change of the 
network topology does not occur during the LDA-HDA 
transition (similar conclusion was also made in Refii^) 
but rather within the HDA phase, between p ~ 2 and 
10 kbar. Concerning the as-prepared HDA', initially be- 



lieved to be the only possible form of HDA, it does not 
seem to have any special importance and represents just 
one particular structure within the HDA megabasin. It 
is not in equilibrium even within the space of amorphous 
structures and its properties are determined by the condi- 
tions of preparation^. As shown above, when prepared 
at T = 80 K, HDA' is rather close to RP at 4.5 kbar. 
The HDA megabasin includes a broad range of struc- 
tures with different local and medium-range order and 
also spans a broad interval of densities. On the high- 
density side, the onset of the VHDA regime marks the 
limit to which the densification can be pushed by adapt- 
ing the network topology, without creating interpenetrat- 
ing networks. The low density limit of HDA stability is 
more difficult to assess with precision. However, as shown 
experimentally by Refs^^^ and also in our simula- 
tions, the HDA region reaches substantially below HDA' 
density. 

From the above it is clear that many forms of HDA ex- 
ist and are metastable at T — 80 K upon decompression 
to p = 0. The important question then is, to what extent 
does this affect the use of the LDA/HDA phenomenology 
as support for the conjecture for the second critical point 
in water. In Refpi& it is claimed that since no unique HDA 
exists, it is difficult to justify the conjecture of a second 
critical point for water. We do not actually think that 
this is necessarily the case. In our picture, the variety of 
HDA ices which exist as metastable forms at p = cor- 
responds to the variety of topologically different HDL at 
different pressures. These various HDA forms cannot in- 
terconvert upon decompression at T = 80 K because this 
would involve rebonding and require overcoming free en- 
ergy barriers. In Ref. 29 it was suggested that VHDA is 
a better candidate for the glassy phase continuous with 
the HDL, rather than HDA'. Since we reserve the use of 
VHDA for the particular high-pressure regime, we can 
say that all forms of the HDA branch of the RP appear 
to be equally good candidates for the glassy phase that is 
the putative continuation of the HDL. Consequently, the 
large density variation of HDA ice at p — might reflect 
the existence of a region of high compressibility of the 
supercooled HDL just below the critical poinlj 2 ^. As sug- 
gested in Refsii^^Si, various forms of LDA ice may also 
exist at T = 80 K and p = 0, but their density variation 
is likely to be much smaller. A sharp transition may well 
exist between the upper limit of LDA and lower limit of 
HDA continua, as suggested by the experiments in Refii^. 
To shed more light on this issue, simulations on larger 
systems combined with finite-size scaling techniques, as 
well as free energy calculations could be helpful. 

In the discussion of the polyamorphism of ice its anal- 
ogy with the properties of silica glasses can be useful, 
although attention should also be paid to important dif- 
ferences. 
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